package statistics;

//                          -*- Mode: Java-Decaf -*- 
// Statistics.java -- 
// Author          : Shiuh Deh Liew
// Created On      : Thu Feb 27 14:55:03 1996
// Last Modified By: Shiuh Deh Liew
// Last Modified On: Sat Apr 13 15:34:08 1996
// RCS: $Id$
//
// Copyright (c) 1995 Shiuh Deh Liew
//
// 

import java.util.*;
import java.lang.*;

/**
 * The Statistics class consists of methods such as Beta1 and Beta0 
 * of regression parameters, average, standard normal distribution,
 * chi-square distribution, t distribution, Simpson's rule, and
 * standard deviation taken from residual sum of square.<BR>
 * @version 1.0 29th February 1996 - Beta1, Beta0, & average
 * @version 1.1 2nd March 1996 added standardDeviation
 * @version 1.2 5th March 1996 added standardNormal
 * @version 1.3 12th March 1996 added tDistribution
 * @version 1.4 12th March 1996 added chiSquareDistribution
 * @version 2.0 6th April 1996 reconstructed Statistics class to meet OOP
 *           standard, and added Statistics instances and regressionStdDev
 * @version 2.1 13th April 1996 added e(residual) for MATH 373
 * @version 2.2 13th April 1996 added correlation and tSignificance
 * @author Shiuh Deh Liew
 */
public class Statistics {

 static private Vector data = new Vector();
 static private Vector dataY = new Vector();


/**
 * Creates an empty Statistics constructor
 * @param none
 * @return Statistics object
 * @exception none
 */
  public Statistics () { }

/**
 * Creates a Statistics constructor that accepts one Vector argument, which
 * contains data.
 * @param Vector
 * @return Statistics object
 * @exception none
 */
  public Statistics (Vector x) { 
    data = x;
}

/**
 * Creates a Statistics constructor that accepts two Vector arguments, which
 * contain data X and Y.
 * @param Vector, Vector
 * @return Statistics object
 * @exception none
 */
  public Statistics (Vector x, Vector y) {  
    data = x;
    dataY = y;
}

  /**
   * Computes Beta1 based on the formular:<P> 
   * Sum X*Y - n*Xavg*Yavg<BR>
   * ----------------------<BR>
   * Sum (X*X) - n*X^2<P>
   * @param int array, int, int array, int
   * @return double 
   */
  public static double beta1 () {
    int sumXY = 0;       //Sum of array X and Y
    int sumXSquare = 0;  //Sum of X^2
    double nXY = 0;      //n * Xavg * Yavg
    double nXSquare = 0; //n * X^2
    int n = data.size(); //data size    
    double x = 0.0;
    double y = 0.0;

    //computes sum of X * Y and sum of X^2
    for (int i=0;i<n;i++) {
      x = Double.valueOf(String.valueOf(data.elementAt(i))).doubleValue();
      y = Double.valueOf(String.valueOf(dataY.elementAt(i))).doubleValue();

      sumXY += x * y;   
      sumXSquare += Math.pow(x,2);
    }
    
    nXY = n * mean() * meanY();
    
    nXSquare = n * Math.pow(mean(), 2);
    
    return (sumXY - nXY)/(sumXSquare - nXSquare);
  }

  /**
   * Computes Beta0 based on the formular:<P>
   * Beta0 = Yavg - (Beta1 * Xavg)<P>
   * @param int array, int, int array, int
   * @return double 
   */
  public static double beta0 () {
    return (meanY() - beta1() * mean());
  }

/**
 * The sum of data
 * @param none
 * @return double
 * @exception
 */
  public static double sum () {
    double sum = 0; //Sum of an array
    int n = data.size(); //data size    

    for (int i=0;i<n;i++) {
      sum += Double.valueOf(String.valueOf(data.elementAt(i))).doubleValue();
    }
    
    return sum;
  }
  
/**
 * The sum of data Y
 * @param none
 * @return double
 * @exception
 */
  public static double sumY () {
    double sum = 0; //Sum of an array
    int n = dataY.size(); //data size    

    for (int i=0;i<n;i++) {
      sum += Double.valueOf(String.valueOf(dataY.elementAt(i))).doubleValue();
    }
    
    return sum;
  }

  /**
   * Calculates the mean of data
   * @param none
   * @return double 
   */
  public static double mean () {
    int n = data.size(); //number of items
    
    return sum()/n;
  }

  /**
   * Calculates the mean of dataY 
   * @param none
   * @return double 
   */
  public static double meanY () {
    int n = dataY.size(); //number of items
    
    return sumY()/n;
  }

  /**
   * Computes Standard Deviation 
   * @param none
   * @return double
   * @exception
   */
  public static double standardDeviation () {
    double SSD = 0.0; //Sum of Squared Deviation
    int n = data.size(); //data size
    double x = 0.0; //temporaty variable to hold data
    double tmpMean = mean();

    //compute the sum of squared deviation
    for (int i=0; i<n; i++) {
      x = Double.valueOf(String.valueOf(data.elementAt(i))).doubleValue();
      SSD += Math.pow((x - tmpMean),2);
    }
    
    return Math.sqrt(SSD/(n - 1));  
  }
  
  /**
   * Computes regression standard deviation
   * @param none
   * @return double
   * @exception
   */
  private static double regressionStdDev () {
    double variance = 0.0;
    double residualSumOfSquare = 0.0;
    double betaZero = beta0();
    double betaOne = beta1();
    int n = data.size(); //data size    
    double x = 0.0;
    double y = 0.0;

    for (int i=0; i<n; i++) {
      x = Double.valueOf(String.valueOf(data.elementAt(i))).doubleValue();
      y = Double.valueOf(String.valueOf(dataY.elementAt(i))).doubleValue();
      residualSumOfSquare += Math.pow((y - betaZero - (betaOne * x)),2);
    }

    variance = residualSumOfSquare/(double)(n - 2); 

    return Math.sqrt(variance);
  }

/**
 * Computes regression interval. It accepts two parameters:<P>
 * t = t distribution<BR>
 * ev = estimated value<P>
 * @param double, double
 * @return double
 * @exception
 */
  public static double regressionInterval (double t, double ev) {
    double sumxxAvgSquare = 0.0;
    double xAvg = mean();
    int n = data.size(); //data size
    double x = 0.0;
    double tmp = 0.0;

    for (int i=0; i<n; i++) {
      x = Double.valueOf(String.valueOf(data.elementAt(i))).doubleValue();
      sumxxAvgSquare += Math.pow((x - xAvg),2);
    }

    tmp= Math.sqrt(1 + 1/(double)n + (Math.pow((ev - xAvg),2)/sumxxAvgSquare));
    return t * regressionStdDev() * tmp;
  }

  /**
   * Standard Normal Distribution. It accepts two parameters:<P>
   * N = number of segments in the numerical integration<BR>
   * x = x-axis<P>
   * @param int, double
   * @return double
   * @exception
   */
  public static double standardNormal (int N, double x) {
    double oldResult = 0.0;
    double E = 0.0001;
    double result = computeStandardNormal(x);
    boolean repeat = true;
    
    if (x > 3.9) return 1.0;
    else if (x < -3.9) return 0.0;

    if (N%2 != 0) {
      System.out.println("N should be an even number.");
      return 0.0;
    }
    
    while (repeat) {
      if ((result - oldResult) <= E) {
    repeat = false;
      } else {
    result = simpsonRule (N, 0, 0, x, 'n');
    oldResult = result;
    N = N*2;
      }
    }

    if (x > 0) {
      result = 0.5 + result;
    } else {
      result = 0.5 - result;
    }

    return result;
  }

  /**
   * Compute Standard Normal Distribution
   * @param double
   * @return double
   * @exception
   */
  private static double computeStandardNormal (double x) {
    return Math.exp(-1 * Math.pow(x,2)/2)/Math.sqrt(2 * Math.PI);
  }


  /**
 * Computes Gamma Function
 * @param int
 * @return double
 * @exception
 */ 
  private static double gammaFunction (int n) {
    double result = 1.0; //The result
    double N = (double)n/2;
    boolean repeat = true;
    
    while (repeat) {
      N = N  - 1;

      if (N <= 0) {
    repeat = false;
      } else {
    result = result * N;
      }
    }
    
    //Gamma(0.5) == Sqrt(Pi)
    if (N%2 != 0) 
      result = result * Math.sqrt(Math.PI);
    
    return result;
  }

/**
 * Computes chi-square Distribution. It accepts three paremeters:<P>
 * N = number of segments in the numerical integration<BR>
 * df = degree of freedom<BR>
 * x = x-axis<P>
 * @param int, int, double
 * @return double
 * @exception
 */
  public static double chiSquareDistribution (int N, int df, double x) {
    double oldResult = 0.0;
    double E = 0.0001;
    double result = computeChiSquareDistribution(df, x);
    boolean repeat = true;
    
    if (N%2 != 0) {
      System.out.println("N should be an even number.");
      return 0.0;
    }
    
    while (repeat) {
      if ((result - oldResult) <= E) {
    repeat = false;
      } else {
    result = simpsonRule (N, df, 0, x, 'X');
    oldResult = result;
    N = N*2;
      }
    }

    return result;
  }
 
  /**
   * Computes chi-square Distribution. 
   * @param int, double
   * @return double
   * @exception
   */ 
  private static double computeChiSquareDistribution (int n, double x) {
    double nominator = Math.pow(x,(((double)n/2) - 1)) * Math.exp(-1 * x/2);
    double denominator = Math.pow(2,((double)n/2)) * gammaFunction(n);
    
    return nominator / denominator;
  }
  
/**
 * t Distribution accepts four parameters, which are:<P>
 * N = number of segments in the integration<BR>
 * df = degree of freedom<BR>
 * x1 and x2 are the integration limits<P>
 * @param int N, int n, double min, double max
 * @return double
 * @exception
 */
  public static double tDistribution (int N, int df, double x1, double x2) {
    double oldResult = 0.0;
    double E = 0.0001;
    double xmax = Math.max(x1,x2);
    double result = computetDistribution(df, xmax);
    boolean repeat = true;
    
    if (N%2 != 0) {
      System.out.println("N should be an even number.");
      return 0.0;
    }
    
    while (repeat) {
      if ((result - oldResult) <= E) {
    repeat = false;
      } else {
    result = simpsonRule (N, df, x1, x2, 't');
    oldResult = result;
    N = N*2;
      }
    }

    if (xmax > 0) {
      result = 0.5 + result;
    } else {
      result = 0.5 - result;
    }

    return result;
  }

  /**
   * Computes t Distribution
   * @param int, double
   * @return double
   * @exception
   */ 
  private static double computetDistribution (int n, double x) {
    double nominator = gammaFunction(n + 1);
    double denominator = Math.sqrt(n * Math.PI) * gammaFunction(n);
    double distribution = Math.pow((1 + Math.pow(x,2)/n),-1*((double)(n + 1)/2));
    
    return nominator * distribution / denominator;
  }

/**
 * Simpson's rule
 * @param int, double, double, double
 * @return double
 * @exception
 */
  private static double simpsonRule (int N, int n, double X1, double X2, char type) {
    double tmpW;
    double tmpTerm = 0.0;
    double result;
    double W;
    double Xlow = Math.min(X1, X2);
    double Xhigh = Math.max(X1, X2);
    double term[] = new double[N];

    tmpW = W = (Xhigh - Xlow)/N;

    switch (type) {
    case 'n':
      term[0] = computeStandardNormal(Xlow) * W/3;
      term[N-1] = computeStandardNormal(Xhigh) * W/3;
      break;
    case 'X':
      term[0] = 0;
      term[N-1] = computeChiSquareDistribution(n, Xhigh) * W/3;
      break;
    case 't':
      term[0] = computetDistribution(n, Xlow) * W/3;
      term[N-1] = computetDistribution(n, Xhigh) * W/3;
      break;
    default:
      System.out.println("Can't find matching distribution type."); 
      return 0.0;
    }
    
    result = term[0] + term[N-1];

    for (int i=1; i<N; i++) {

      switch (type) {
      case 'n':
    tmpTerm = computeStandardNormal(tmpW);
    break;
      case 'X':
    tmpTerm = computeChiSquareDistribution(n, tmpW);
    break;
      case 't':
    tmpTerm = computetDistribution(n, tmpW);
    break;
      }

      if (i%2 == 0) {
    term[i] = 2 * tmpTerm * W/3;
      } else {
    term[i] = 4 * tmpTerm * W/3;
      }
      tmpW += W;
      result += term[i];
    }

    return result;
  }

  /**
   * y hat = beta0 - beta1*x
   * @param double
   * @return double
   * @exception
   */
  public static double yHat (double x) {
    return beta0() + (beta1() * x);
  }
  
  /**
   * e - residual (y - y hat) 
   * @param none
   * @return double
   * @exception
   */
  public static void e () {
    int n = dataY.size();
    double e = 0.0;
    double x = 0.0;
    double y = 0.0;
    
    for (int i=0;i<n;i++) {
      x = Double.valueOf(String.valueOf(data.elementAt(i))).doubleValue();
      y = Double.valueOf(String.valueOf(dataY.elementAt(i))).doubleValue();
      e = y - yHat(x);
      System.out.println("x: " + x + " y: " + y + " e= " + e);
    }
    
  }

/**
 * SSE - Error Sum of Square 
 * @param
 * @return
 * @exception
 */
  public static double SSE () {
    int n = dataY.size();
    double sse = 0.0; //error sum of square
    double x = 0.0;
    double y = 0.0;
    
    for (int i=0; i<n; i++) {
      x = Double.valueOf(String.valueOf(data.elementAt(i))).doubleValue();
      y = Double.valueOf(String.valueOf(dataY.elementAt(i))).doubleValue();
      sse += Math.pow((y - yHat(x)),2);
    }
    
    return sse;
  }

/**
 * standardError - Standard error of the estimate
 * @param
 * @return
 * @exception
 */
  public static double standardError () {
    int n = dataY.size();
    
    return Math.sqrt(SSE()/(double)(n - 2));
  }

/**
 * correlation
 * @param
 * @return
 * @exception
 */
  public static double correlation () {
    int n = data.size();
    double x = 0.0;
    double y = 0.0;
    double sumX = 0.0;  //sum of x
    double sumY = 0.0;  //sum of y
    double sumXY = 0.0; //sum of x*y
    double sumXSquare = 0.0;  //sum of x square
    double sumYSquare = 0.0;  //sum of y square
    double numerator = 0.0;   //numerator of correlation
    double denomerator = 0.0; //denomerator of correlation

    //computes sum of ...
    //sum of x and y are computered here to improve performance; eventhough,
    //there are methods for sum of x and y.
    for (int i=0;i<n;i++) {
      x = Double.valueOf(String.valueOf(data.elementAt(i))).doubleValue();
      y = Double.valueOf(String.valueOf(dataY.elementAt(i))).doubleValue();
      
      sumX += x;
      sumY += y;
      sumXY += x * y;   
      sumXSquare += x * x;
      sumYSquare += y * y;
    }    

    numerator = n * sumXY - sumX * sumY;
    
    denomerator = Math.sqrt((n*sumXSquare - Math.pow(sumX,2)) *
      (n*sumYSquare - Math.pow(sumY,2)));

    return numerator/denomerator;
  }

  /**
   * tSignificance
   * @param none
   * @return double
   * @exception
   */
  public static double tSignificance () {
    int n = data.size();    
    double r = correlation();

    return (Math.abs(r) * Math.sqrt(n - 2))/(Math.sqrt(1 - Math.pow(r,2))); 
  }

}
